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I. INTRODUCTION 

t-H \ 

The study of the density matrices of identical particles (bosons or fermions) moving freely in a box jlj is generalized 
in this paper to the case of identical particles in a parabolic confinement potential with either harmonic interactions 
between the particles, or with an anisotropy induced by a homogeneous magnetic field on top of the parabolic 
, confinement. This model, giving rise to repetitive Gaussian integrals, allows to derive an analytical expression for 
the generating function of the partition function. For an ideal gas of non-interacting particles in a parabolic well, 
this generating function coincides with the grand-canonical partition function. With interactions, the calculation 
of this generating function circumvents the constraints on the summation over the cycles of the permutation group. 
Moreover, it allows one to calculate the canonical partition function recursively for the system with harmonic two-body 
interactions. The theory is developed both for fermions and for bosons. In view of the recent interest in Bosc-Einstein 
condensation in a trap more attention has been payed to the boson case in the discussion of the results. The 

model system discussed here has already been studied in the context of quantum dots with operator techniques, 
and the eigenvalues and eigenstates were calculated including the effect of harmonic two-body interactions and in the 
presence of a magnetic field ||. However, to the best of our knowledge neither the boson case nor the thermodynamics 
O ■ seems to have been analyzed previously. It should also be mentioned that the idea of first expanding the Hilbert space 
to the configuration space and then projecting onto the appropriate subspace by group theoretical means has been 
used recently HQ in the context of quantum dots to study the ground state correlations for fermions and bosons. 

Another motivation to perform the present analytical calculations lies in our path integral formulation of the density 
matrices for N particles []s|-|lO|]. Indeed, particles in a parabolic potential are a favorite testing ground for the path 
integral method JlT|-p^|. It should be noted that for our formulation the existence of a positive measure over a 
well defined domain in the R 3N configuration space is essential in view of any algorithmic approach to the problem. 
In the present paper, which is in essence analytical, integrations over the configuration space are performed. The 
reason is that the extension of the state space to the configuration space makes the Gaussian integrals tractable. 
The permutation symmetry leads to summations over the cycles which are performed using the generating function 
technique, which is one of the main results of the present paper. 

The model of N identical particles in a parabolic well, in the presence of a magnetic field and with harmonic 
repulsive or attractive two-body interactions has its intrinsic value, since it constitutes an exactly soluble idealization 
of atoms in a magnetic trap. It should be stressed that the association of identical particles with each three oscillator 
degrees of freedom makes the model 3D. Without Bose-Einstein or Fermi-Dirac statistics, i.e. for "distinguishable" 
particles, the model is equivalent with 3N one dimensional oscillators, because each degree of freedom decouples in 
such a way that there is no difference in statistical behavior between 3N ID oscillators and N 3D oscillators [Q. 

This paper is organized as follows: the calculation technique is explained in the next section. In the subsequent 
section we repeat the same calculation for the model with a homogeneous magnetic field. In section IV the results 
for ID bosons and fermions are given. In section V bosons in 3D are analyzed in some detail, and in the last section 
the conclusions are given. 
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II. HARMONICALLY INTERACTING IDENTICAL PARTICLES IN A PARABOLIC WELL. 



In this section we calculate the partition function of N identical particles with the following Lagrangian including 
one-body and two-body potentials 



L = \Y,*}- V i- V * * = T^> 2; ^-^E^-rO 2 . (2.1) 

3=1 3=1 3,1=1 

(Atomic units are used.) The potentials can be rewritten in terms of the center-of-mass coordinate R and coordinates 



Uj describing the coordinates of the particles measured from the center of mass 

N 



3=1 

from which 



R= ^E r ^ u^rj.-R, (2.2) 



1 N 
V 1 + V 2 = V CM + V; V CM = -NtfR 2 : V = w 2 ^ V, ( 2 - 3 ) 



2 

with 



3=1 



w = ^n 2 - Ncu 2 . (2.4) 

The requirement that w has to be positive expresses the stability condition that the confining potential has to be 
strong enough to overcome the repulsion between the particles. If an harmonic interparticle attraction is considered, 
the eigenfrequency w would become w — \/£! 2 + Nlu 2 , and no stability condition has to be imposed on the confining 
potential. Notice that these transformations do neither diagonalize the Lagrangian nor the Hamiltonian, because the 
coordinates Uj are not independent of the center-of-mass coordinate. 

Since the system consists in each direction of 1 degree of freedom with frequency fl and (N — 1) degrees of freedom 
with frequency w, the propagator 

K D (r'{ ■ ■ ■ r%, flri ■ • -r^O) = <r'/ • • • r% |e"^| r[ ■ ■ ■ r' N ) D (2.5) 

for distinguishable particles (indicated by the subscript D for 3 dimensions and d in 1 dimension) can be calculated 
from the action expressed in the imaginary time variable, and it is of course a product of the propagators Kd per 
component: 

K D {r'{ ■ ■ ■ r&, I3\r[ ■ ■ ■ r' N , 0) = K d (x", (5\x' , 0) K d (y", f3\y' , 0) K d (z", (3\z', 0) , (2.6) 

where the column vector x contains the ^-components of the particles, i.e. x T — (xi, • • • , xm) and similarly for y and 
z. Knowing the propagator K (x", (3\x' , 0) CT of a single harmonic oscillator with frequency w 



(^x 2 p + Xq^J cosh w(5 — 2x/3Xq 



K ^M m = ] /^— ^exp j--^ ^ } . 

one finds for the 1-dimensional propagator K d of the N distinguishable oscillators in the interacting system: 

K (^X",/3\y/NX',Q) n 
K d (x", (3\x', 0) - -) 7 =— i - — J] K (4, P\x' v Q) w , (2.8) 

3=1 



K (yNX",f3\y/NX',0 



where the factor v N in v NX" accounts for the mass N (in atomic units) of the center. The denominator in ( |2.8| ) 
compensates for the fact that (N — 1) instead of N degrees of freedom of frequency w are availabl e. T he 3-d ime nsional 
propagator K D (^1]) for N distinguishable oscillators of the interacting system is, according to (2.6) and (|2.8D, given 
by 
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JY 



K (VNR",[3\VN~R.',o) 
K D (f",/3|r', 0) = -) — - — ^ [] AT (r?, (3\r' J7 0) w , (2.9) 



K \yNBJ' , (3\\f N~R! ,Q 

where r denotes a point in the configuration space R 3N , i.e. f T = ( (xi, y%,Zi) , • • • , (aiiv, yjv, zjv) j . The symmetrized 
density matrix if/ for 3D identical particles (indicated by the subscript /) can be obtained by using the following 
projection, with P denoting the permutation matrix: 

Kj (r", /3|r', 0) = ^ £ (Pf, /3|f 0) , (2.11) 

p 

where £ = +1 for bosons and £ = — 1 for fcrmions. It should be emphasized that P acts on the particle indices, not on 
the components of r separately. The partition function is then readily obtained by integrating over the configuration 
space 

Zj = f dfKj (r, /3|r, 0) = f dfj- £ ^K D (Pf, /3|r, 0) . (2.12) 

The remaining part of this section will be devoted to the explicit evaluation of this integral for the partition function. 
The integration proceeds in 3 stages: the first stage deals with the center-of-mass treatment, the second one concerns 
the cyclic decomposition, and in the third step the summation over the cycles will be performed. 



A. The center of mass 



The center-of-mass coordinate R does not only depend on the coordinates of all the particles, but it also has its 
own propagator. Therefore substituting R by its expression in terms of the particle positions and then performing 
the integration seems not to be the most adequate way to deal with the integration over the configuration space. 
Instead, the following identity is used for the formal treatment of R as an independent coordinate, at the expense of 
additional integrations: 

/ df f (^>jf E =J dR l dff ( F ' R ) 6 jf E r ^ ■ ( 2 - 13 ) 



Fourier transformation of the <5-function then leads to 

N 




dk 

(2nf 



dR, j 7^3e lk R / dff (r, R) e~ ik ' ? , (2.14) 



where k T = -^((1,1,1), (1,1,1)) is a 3 AT dimensional row vector. Applying this transformation to the partition 

function Zi and rearranging the factors one obtains 



r dv k(Vnr,/3\Vn-r,o) , 1 n ^ 

dR / -^e ikR — ^ / ^E^n^f^'^i' ) ( 2 - 15 ) 



(2tt) 3 k (VNR, f3\ VNTi, 



P 3 = 1 



This transformation makes R independent of the particle positions relative to the center of mass. The real de- 
pendence on the relative positions is reintroduced by the Fourier transform. It should be noted that the explicit 
dependence of the propagator (2.9) on R, and the presence of the factor e~ lk ' rj / N , are consequences of the two-body 
interactions. 

The next step is to rewrite the sum over the permutations as a sum over all possible cycles. This will be done in 
the next subsection. An excellent example of such a decomposition into cycles has been given by Feynman Q| for a 
system of non-interacting particles in a box. 
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B. Cyclic decomposition 



A permutation can be broken up into cycles. Suppose that a particular permutation contains Mi cycles of length I. 
The positive integers Mi and I then have to satisfy the constraint 

J2?Mi = N. (2.16) 

i 

Furthermore, the number M (Mi, • • • Mat) of cyclic decompositions with Mi cycles of length 1, • • • , Me cycles of length 
I, ■ ■ ■ is known to be 

^•••^) = n^v (2 - 17) 

A cycle of length I will be obtained from [i — 1) permutations. Therefore the sign factor £ p can be decomposed as 

£P=JJ^-i)M f _ (218) 

Combining these result originating from the permutation symmetry one obtains 



Zl = dR -^V k '*_) (o. y TT S-_ (/Q (k)) Mf , (2.19) 



N 

K t (k)= /" dr,+i Jdxf- J dr 1 5{vi +1 ~v 1 )\[K{r ]+u l3\v 3 ^) w e- i 

The 5-function expresses that the decomposition is cyclic. It is obvious that 

K t (k) = /C< 11?) (k x ) K { e 1D) (k y ) 4 1D) (k z ) , 
which allows to analyze Ki (k) from its 1-dimcnsional constituents: 

K-i 1D) (k x ) = dxi +1 dxf- dxi5(x i+ i-x 1 )Y\_K(x j+ i,f3\x 3 ,0) w e 

j=i 



k-rj/N 



-ik x Xj /N 



(2.20) 



(2.21) 



(2.22) 



Using the semigroup property |l5[] of the harmonic oscillator propagator K (iCj+i, /3 1 2;^ , 0) , all integrations but one 
can be performed 



K, ( l D) {k x ) = / dxK(x,^\x,O) w e-Io ^f^Mr) 



(2.23) 



where 



3=0 



The integral ( 2. 23] ) is the propagator K w j of a driven harmonic oscillator with the Lagrangian 



L w j x = l;x 2 - T^ 2 ^ 2 + fx (t) x. 



(2.24) 



(2.25) 



studied in f pd|JT^ ]. It should be noted again that without two-body interactions the driving force ( [2.24 ) is lacking. 
Taking over the result from |16| and integrating over the configuration space one obtains 
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Zw,f x (P) = J dxK w j x (x,f3\x,0) 



1 



2sinh ±/3w 



exp 




P d JAr)fA<y) cosh^f-|r-a|j W j 



2(0 



After straightforward algebra one obtains for the ID function K^ D ^ (k x ) 



(ID) 



and for its 3D extension: 



K, 



K-e ( k 



(k x ) 



2 sinh \k 



1 

2 sinh \if3w 



■ exp 



exp 



I k 2 x l + e 
' 4N 2 w 1 - e~P w 



sinh i/3u> 



-/3iu 



£ fc 2 1 



-/3u> 



AN 2 w 1 



Using ( [2.16 ) one then is left with a sixfold integral for the partition function 



k(Vnr,/3\Vnr 



l^M t --M N iU M £ \e M t \ 2 sinh | 



Jf^ e -*-P ^ 4jV a, l-e-f™ y ■ 
\ 3Mi 



■M N ll£ M f !« M * V 2sinh 3 £ / 3 "', 

Both the integrations over k and R are Gaussian, leading to the following series for Zj: 



Zj = 



( sinh i/3w 
Vsinhi/30 



Z z (iV); Zj(JV)= ^ [] 

Mi — M w f 



(*-l)M< 



1 _ p-tftw 



3M ( 



(2.26) 



(2.27) 



(2.28) 



(2.29) 



(2.30) 



Without two-body interactions (w = fi), Zj (AT) is the partition function of a set of identical oscillators. The 
partition function Zj only differs from it by a center-of-mass correct ion a nd the actual values of w. 

The remaining summation over the cycles involves the constraint (2.16), which however can be removed by the use 
of the generating function technique, which will be considered in the next subsection. 



C. The generating function 

Concentrating on the explicit dependence of Zj (N) on N (with w considered as a parameter), one can construct 
the following generating function 

oo 

S(u) = %i(N)u N , (2.31) 

with Zj (0) = 1 by definition. The partition function Z/ (N) can then be obtained by taking the appropriate derivatives 
of 5 (u) with respect to u, assuming that the series for 2 (u) is convergent near u = 0: 



u— U 

The summation over the number of cycles with length i is now unrestricted, and can easily be performed: 



(2.32) 



This scries can be rewritten into the more familiar form 



Hz («) = exp (-£ £ i [y + 1) [y + 2) In (l - £ ue -* u (*+''))^ . (2.34) 

It should be noted in view of the remarks at the end of the preceding section concerning Zj (N), that 2/ (u) is also 
the generating function of a model without two-body interactions. In that case w equals fi, and Sj (u) coincides with 
the well known [jl7|-|l9| grand canonical partition function of a set of identical particles in a parabolic well. 
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D. Recurrence relations for the partition function 



Starting from the expression for 3/ (u) derived in the previous subsection for the interacting model, a recursion 
relation can be obtained for Zj (N). Introducing: 

b = e~^ w (2.35) 

for brevity in the notations, we observe that 

d °° 1 b^"^ 11 

-Zt (u) = («) X - („ + l) („ + 2) _ ^ (2.36) 



i/=0 



Considering next Zj (iV) = m d^jv-i SI"( u )L-o ! ^ e P r °duct rule and an elementary binomial expansion can be 
used to find 

1 N ^ ( hh(N-m) \ 3 

m=0 \ / 

The corresponding 1-dimcnsional version of this recurrence relation (indicated with the subscript i to distinguish 
it from the 3D case with capital subscript) becomes 

^ (^) = jvE f^Tzw^ 1 * {m) ' (2 - 38) 

m— 

leading to the following partition functions in closed form for one-dimensional bosons and one-dimensional fermions 

bi N bi N ' 
Z b = — » ; Z/ = —f, . (2.39) 

nf=i(i-^) nf =1 (i-^) 

It is easy to check that these partition functions are the solution of the recurrence relation for Zj (N) with £ = 1 for 
bosons and £ = — 1 for fermions. However we did not find a systematic method to obtain analytical solutions of this 
type of recurrence relations. E.g. for the 3D case we had to rely on numerical schemes, as will be discussed below. 
But at this stage, it is worthwhile first to consider the presence of an homogeneous magnetic field, as the origin of 
anisotropy in our model of TV identical oscillators. 

III. N IDENTICAL OSCILLATORS IN A MAGNETIC FIELD 

The Lagrangian of N particles in a confining parabolic potential in the presence of a magnetic field is 

N N 

= 2 L ft - tocXjVj) 2 - if 2 £ r|, (3.1) 

J=l 3=1 

where ui c is the cyclotron frequency. For this model, the calculations of the preceding section can in essence be 
repeated. First the propagator for distinguishable particles is calculated. The next step will be the projection on the 
irreducible representation of the permutation group, and performing the cyclic decomposition. Then the generating 
function is introduced to circumvent the constraints on the partition in cycles, and finally the summation over the 
cycles is performed. 

The fact that the energy spectrum and the wavefunction can be calculated when harmonic interparticle interactions 
are included S indicates that the propagator and the partition function for the model in a magnetic field with 
two-body interaction can be obtained using our methods. The calculation technique is very demanding, therefore it 



seemed appropriate to illustrate the method at the hand of the simple model (3.1) 



G 



A. The propagator for distinguishable particles 



For distinguishable particles, the many-particle propagator is a product of one-particle propagators. The one- 
particle pr opag ator K^} (r, /3|r') of this model can be calculated by stochastic techniques pfj] or by path integral 
techniques (ill. The evaluation is somewhat lengthy but straightforward, resulting eventually in 



n((z 2 + (z') 2 ) cosh/3n-2zz') 
2tt sinh /3£2 2tt sinh (3s eX P i 2 sinh /3S1 

{„ (x 2 +y 2 + (x') 2 + (y') 2 ) cosh /3s— 2(a:a: / +yy / ) cosh 
— o " ■ i, a 

2 sinh ps 

x exp |-i (±w c (an/ - x'y') - s^j^ (y'z - yx'yj } , 



(3.2) 



with the eigenfrequency s given by 



fi2 + ^ c 2 . 



(3.3) 



The partition function corresponding to this single-particle propagator is obtained by integrating over the configuration 
space: 



Z&> 09) = j drKi 1 } (r,/3|r) = 1 8 sinh sinh /3 (| + ^) sinh/3 (| - ^) 
which coincides with the partition function of a 3D harmonic oscillator if oj c = 0. 



(3.4) 



B. The partition function and generating function for identical particles 

For N identical particles, the propagator becomes 



p 3=1 

with the corresponding partition function given by 



1 N 

K I<Uc (r, P\T>) = M E^II K % (t Pr h > ^i) > ( 3 - 5 ) 



Zj,o, e (AT) = / drA7, Wc (r, f3\f) . (3.6) 

Using the cyclic decomposition and the semigroup property as before, followed by the relaxation of the constraint 
(2.16) on the number of cycles, one finds for the generating function 



\C=1 



e-i 



s Wc («) = ex P rr m u < . (3.7) 



Because in this model with a magnetic field no two-body interactions were taken into account, this expression equals 
the grand-canonical partition function of the model if u = is interpreted as the fugacity, with fx the chemical 
potential. 

We have thus obtained an expression for the grand-canonical partition function in the absence of two-body inter- 
actions for N identical harmonic oscillators, both with and without a magnetic field. For a system with harmonic 
two-body interactions, repulsive or attractive, expressions for the canonical partition function are obtained (in the 
absence of a magnetic field) . The corresponding thermodynamic properties (grand-canonical ground state occupancy, 
canonical and grand-canonical specific heat and average square radius) of the confined system will be calculated in 
the next section. 
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IV. IDENTICAL PARTICLES IN ID 



In the preceding sections we have obtained the grand canonical partition func tion of identical harmonic oscillators 



in a confining para bolic potential well without a magnetic field [equation (2.34)] and in the presence of a magnetic 
field [equation J3.7| )1. For the same system the canoni cal pa rtition function is obtained for an attractive as well as for 
a repulsive harmonic two-body interaction [equation ( [2.30 )]. The expressions are given for bosons and for fermions. 



In this section we will analyze the thcrmodynamical properties of the confined system in the case of motion in ID. 
It should be noted that for the non-interacting case in ID, we recover the results of Ref. To the best of our 

knowledge, the interacting case has not been analyzed up to now for identical particles. For distinguishable particles 
it has been studied before p2|| . 

In ID, the explicit expressions for Z largely facilitate the analysis: the free energy F = — i huZ, the internal energy 

U = -Jg- and the specific heat C = ^ can be readily obtained. The results are summarized in Table I. The subscript 
CM refers to the center-of-mass contribution. 

The center-of-mass correction Ccm to the specific heat -which is an excess specific heat and therefore may be 
negative- is shown in Fig. 1. In Fig. 2 we plot the contribution AC to the specific heat, due to the internal degrees 
of freedom. Note that AC is the same for ID fermions and bosons; bosonization implies this result |23[| . 

The spatial extension of the cloud of interacting identical particles can be obtained by applying the Hellman- 
Feynman theorem to the free energy: 

^-s*(--"5>s)-7<**>- 

The difference \X/f — \ x /f, between the mean square length of the ID fermion and boson system turns out to be: 

<* 2 >/-<A~- (4-2) 

It is proportional to the number of particles, and inversely proportional to the frequency w of the internal degrees of 
freedom. 



V. BOSONS IN 3D 



The calculation of the thermodynamic properties or the mean square radius of the 3D cloud of identical particles 
is substantially complicated by the lack of an explicit expression in closed form for the partition function. Numerical 
methods have therefore to be used. We concentrate this analysis on the boson case. For fermions a different scheme 
will have to be developed and will be published later. 



A. In the absence of a magnetic field 

The recurrence relation ( 2.37| ) is not directly accessible for numerical computation, as can easily be seen by evaluating 
the expected dominant factor b 3N / 2 / Y[jLi (l — &) 3 for bosons. For a relatively low temperature and a moderate 
number of particles, say b = 0.75 and N = 1000, this factor is as small as 1.0402 x 10~ 182 . We therefore isolate this 
factor using the following scaling for the partition function 

Z B (N)=a N 2 (5.1) 

n,=i(i-^) 



and rewrite the recurrence relation (2.37) in terms of the activity pn (see e.g. |22j) defined as 



N 

ctn = Pn<Jn-\ => cr 7v = <t rj Pj (5.2) 

3=0 



where do = Po — 1 have been introduced for convenience. The recurrence relation (2.37) then becomes after some 
manipulations 
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1 /l-^N 3 / l_ 6 \3 iV-1 



1+ e n ^ • m 



AT V 1 - 6 7 I ^ V 1 - 6^-™ 7 11 p,- 

The corresponding recurrence relation for the internal energy of the internal degrees of freedom becomes 

3 1+6 , Vb(N-I) 

U b (JV) 1 1 (1- b N \ a I 2I=E + — mL 



hw N 



1 fl-b N \ 3 2—+ hw \ 

^{—) ^ +E »-((i (w - m){ ^ + !!^)(^) s n^ +I i!#)J- (M 



The temperature scale used to express b — e w / kT is 



N \ 1/3 kT T 



t= ^W)J — = ^ where C (3) = 1-2021. (5.5) 



The routine used to calculate the values p± ■ ■ ■ pn and Ub (1) • • - Us (N) from these recurrence relations (5.3) and 



(5.4) is displayed in Table II. 

The values of p\ ■ ■ ■ pn for N = 100 are shown in Fig. 3 for the three temperatures T/T c = 0.5, 1 and 2. For 
T = 2T C it turns out that pj is about 8 times larger than Pj-\ for j approaching N, which makes it extremely difficult 
to deal numerically with the values of the partition function rather than with the proportionality factors. 

Once pn is known, the internal energy Us (N) is readily obtained from (^4|). The results for N = 10, 100 and 1000 
are shown in Fig. 4. The specific heat in the canonical ensemble can also be calculated from these parameters. This 
is shown in Fig. 5, clearly illustrating the effect of condensation for a finite number of particles. Anticipating the next 
subsection, it should be noted that the specific heat in the canonical ensemble for N particles and without repulsive 
two-body interactions is identical to the specific heat in the grand canonical ensemble with the average number of 
particles given by the same value of N. The reason why we compare both ensembles without repulsive interactions is 
the requirement that in the grand canonical ensemble the system should be stable for any large number of particles, 
which is not the case in the Gaussian model wit h re pulsion, because the confining potential can only accommodate 



a finite number of particles as a consequence of (2.4). Another interesting consequence of the repulsive interactions 



follows from the de pend ence of the condensation temperature T c on the number of particles, which obeys the following 



scaling law taking (2A) into account: 



C(3)y n V fi J V n 2 \ n 

For the case of attractive interactions and no confinement potential, T c is proportional to iV 4 / 3 . The condensation 
temperature for both cases is plotted in Fig. 6. For the Gaussian model, the case of harmonic attraction does not 
pose any problem because in this model there is no sign of an "extra" collapse due to the nature of the interaction. 
The only consequence seems to be that the condensation occurs at a much higher temperature than would be the 
case without two-body interactions or with a repulsive harmonic interaction. 

B. In the presence of a magnetic field 

The thermodynamic properties in the grand canonical ensemble of an ideal Bose gas in a parabolic well and in the 



presence of a magnetic field can be obtained directly from ( |3~4| ) and (3.7). Substituting u by the fugacity e^ M , the 
Gibbs free energy G Wc for the boson case (£ = +1) becomes: 



U / (1 - e-*°) (l - e -M*+^)) (l - „-*(-*«.)) ' (5J) 
After a power series expansion of the denominators, the summation over the cycle lengths t can be performed: 



G ^ = \ £ hx(l-e^e- 0n (i +j )e-^ 1+k + l ^e-^ k - 1 ^, (5.8) 
as expected from a Bose-Einstein distribution with energy levels 
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e j,k,i = ^Cl + s + jQ + k ys + ^uicj + l (s - j . ( r,::) i 

dG 

The average number of particles N = — is given by 

N = ) nj k,i\ nj,k,i = T-77 , a ; ; — (5.10) 

.fji Q J 1 _ e /3Me"^(2+j) e -/3 S (l+fc+0 e -5^<=( fe -0 

and the fugacity can be eliminated in favor of the ground state occupancy no: 



n 



^ = ae 0(tn +s ) with a s (5J1) 
1 _ ePi* e -hPn e -ps n + l 

Restoring the cyclic summation for reasons of numerical convergence, one obtains for the average number of particles 

°° „4 



N 



£fL_ ; Dl = (i _ e -^ n )(l - e -^(*+i-)) (l - e -^-K> ) , (5.12) 



which can be cast in the form of the following numerically tractable series 



» = i^ + £"'GH- (5 - 13) 

For given N, standard numerical techniques can be used to determine a € [0, 1] , and hence the ground state occupancy 
no, which is shown in Fig. 7 for lu c — and in Fig. 8 for u> c /Vt = 5 as a function of T/T c for several values of N. It 
turns out that the parabolic well is more important than the anisotropy due to the presence of the magnetic field. 
The anisotropy only moderately influences the temperature dependence of the ground state occupancy. 
The internal energy U Uc = — /aN becomes 

^ = g ^ (°T3^ + (- + 2^ j TTTM^) + I 5 " 2 WC J x-e-^-^) J (5 ' 14) 

and its numerical evaluation once a is determined presents no numerical difficulties. The resulting specific heat 
is shown in Fig. 9 for w c /Q = 5. Comparison with Fig. 5 reveals that the anisotropy due to the magnetic field 
essentially broadens the peak in the specific heat near the condensation temperature, but does not substantially alter 
the structure. 



VI. DISCUSSION AND CONCLUSION 



In this paper we applied the method of symmetrical density matrices, developed by Feynman for a system of non- 
interacting particles in a box, to a system of harmonically interacting particles in a confining parabolic potential. The 
interaction could be taken into account, due to the Gaussian nature of the propagators, allowing integration over the 
configuration space. The symmetrization resulting from the projection of the propagators for distinguishable particles 
on the appropriate representation of the permutation group gives rise to a series which could be summed using the 
generating function technique. Without these generating functions the calculation has to be restricted to a limited 
number of particles . Using them not only the grand canonical partition function 5 (u) could be obtained in 

the parameter range of the model where S (u) is well defined, but also the canonical partition functions Z (TV) for 
a given number N of identical particles could be obtained as a recursion of partition functions of a smaller number 
of particles for the interacting system. The recurrence relation for the activity pjv, i.e. the proportionality factor 
between Z (TV) and Z (N — 1) , allows for an accurate numerical treatment of the thermodynamical quantities of the 
model, such as the internal energy and the specific heat. 

Also the thermodynamic properties of the same model in the presence of a homogeneous magnetic field could be 
investigated along these lines. A detailed analysis of the ideal gas in a confining parabolic potential with anisotropy 
due to a magnetic field was presented; special attention was payed to the relation between the number of condensed 
atoms, the magnetic field, the strength of the confinement potential and the total number of particles. The relationship 
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between the parameters of our model and the characteristics of atomic traps [£]-§} lies beyond the scope of the present 
paper. 

It should be mentioned that in the absence of two-body interactions, the generating function can be identified as 
the grand canonical partition function. In that case the specific heat as obtained from the grand canonical ensemble 
for an average number of (N) particles calculated using the chemical potential equals the specific heat obtained from 
the canonical ensemble with the number of particles iV given by (N) . For a more elaborate discussion we refer to |l7| . 

It should also be mentioned that in this the model it is assumed that the spin degrees of freedom are fixed. This 
simplifying assumption is imposed by the symmetrization method, which becomes more involved if the spin degrees 
of freedom depend on the configuration of the particles. 

In summary, the partition function of a general Gaussian model for bosons and fermions, with or without a magnetic 
field, has been calculated analytically, and the thermodynamical properties of this model have been studied. Particular 
attention has been given to the Bose-Einstein condensation in the presence of two-body interactions, repulsive and 
attractive, and in the presence of a magnetic field. 
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Tables 



Center-of-mass contribution 



/ (fermions) 



i = b (bosons) 



Z = ZcM^i 



iinhm/3/2 
iinh {2/3/2 



/2 



nCT 



nf = i('-- w ) 



F = Fcm + Wi 



2 ^ a m i- 



\N 2 w + AF 



±Nw + AF 



U = Ucm + V, 



\N 2 w + AU 



iNw + AU 



AU - 



C = Ccm + Cj 



(1- 



-fl/3 



) 2 



AC 



AC 



AC = fc/3 2 w 2 J] 



' =1 (i- 



TABLE I. Partition function, free energy, internal energy and specific heat of ID identical particles in a parabolic confining 
potential with harmonic two-body interactions. 



subroutine recur (MAX , b , RHD , U) 
dimension RHD (0: MAX) ,U(0:NMAX) 
RH0(1)=1 

U(l)=1.5*(l+b)/(l-b) 
if (NMAX.le.l) return 
do N=1,NMAX 
RH0(N)=0 
do m=0,N-2 

RHD (N)=(l-b** (m+1 ) ) **3/RH0 (m+1) * 
k ( RH0(N)+( (l-b)/(l-b**(N-m)) )**3 ) 

enddo 

RH0(N)=(1+RH0(N))/N * ( (l-b**N) / (1-b) ) **3 
U(N)=0 

do m=0,N-2 

DU=1.5*(N-m)*(l+b**(N-m))/(l-b**(N-m)) + U(m) 

U(N)=(l-b*>Km+l))**3/RH0(m+l)* 
* ( U(N)+DU*( (l-b)/(l-b**(N-m)) )**3 ) 

enddo 

U(N)=l./N/RH0(N)*((l-b**N)/(l-b))**3* 
•= (U(N)+1.5*(l+b)/(l-b)+U(N-l)) 
enddo 
end 



TABLE II. Subroutine to calculate the values pi • • • pjv and Ub (1) • ■ - Us (iV) from the recurrence relations (5.3) and (5.4). 
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Figure captions 

Fig. 1: Center-of-mass correction Ccm to the specific heat of ID fermions and bosons (see table I) in units of the 
Boltzmann constant k as a function of the reduced temperature tq = kT/Ct, for various ratios of w/Q. 

Fig. 2: Specific heat contribution AC/Nk per particle for ID fermions and bosons (see table I) in units of the 
Boltzmann constant k as a function of t/N with r the reduced temperature r = kT/w, for several values of N. 

Fig. 3: Values of pi ■ ■ ■ pn for N = 100 from the recurrence relation ( |5.3| ) for the three temperatures T/T c = 0.5, 1 
and 2. 

Fig. 4: Internal energy u/N = V B (N) /Nhw for N = 10, 100, 1000 as a function of T/T c . 

Fig. 5: Specific heat Cs/Nk per particle in units of the Boltzmann constant k for N — 10, 100, 1000 as a function 
of T/T c . 

Fig. 6: Critical temperature T c as a function of the number of particles N for a repulsive (full line) and an attractive 
(dashed line) two-body interaction. 

Fig. 7: Relative ground state occupancy no/N as a function of t/N 1 / 3 with r the reduced temperature r = kT/w 
for uj c = and for several values of N. 

Fig. 8: Relative ground state occupancy tlq/N as a function of t/N 1 / 3 with r the reduced temperature r = kT/w 
for o; c /f2 = 5 and for several values of iV. 

Fig. 9: Specific heat Cs/JVfc per particle in units of the Boltzmann constant A: as a function of t/N 1 / 3 with r the 
reduced temperature r = kT/w for u c /rt = 5 and for several values of TV. 
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